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In this letter we discuss use of Granger causality to the analyze systems of coupled circular 
variables, by modifying a recently proposed method for multivariate analysis of causality. We show 
, the application of the proposed approach on several Kuramoto systems, in particular one living 

• on networks built by preferential attachment and a model for the transition from deeply to lightly 

anaesthetized states. Granger causalities describe the flow of information among variables. 
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Countless physical systems are effectively represented using circular variables, such as phases or orientations [l|. 
£C) . Coupled oscillators systems describe, e.g., arrays of Josephson junctions, chemical reaction diffusion systems, circa- 
t-H ' dian oscillations, brain activity, and many others 2]. A problem of particular interest is to assess the interaction 
between sub-systems, each described by a phase variable, a task which has been tackled using the ideas of generalized 
' (— i '■ synchronization, phase synchronization and phase dynamics modeling [3]. A great deal of attention has been recently 
paid to the interplay between the properties of the coupling network and synchronization of oscillators [i[ . 

Granger [5| proposed a major approach to analyze causality between two time series: if the prediction error of 
the first time series is reduced by exploiting the knowledge of the second one, then the second time series is said 
to have a causal influence on the first one. Initially developed for econometric applications, Granger causality has 



to nave a causal mnucncc on tnc nrst one. initially developed lor econometric applications, Granger causality nas 
gained popularity also among physicists (see, e.g., [a, 0, H, HI ) ■ Being closely related to transfer entropy [To!) . Granger 
causality is connected to the amount of information being transferred from one time series to the other. A novel 
approach for multivariate Granger causality has been recently proposed in [TlT |: the problem of false-causalities is 
\ addressed by a selection strategy of the eigenvectors of a reduced Gram matrix whose range represents the additional 
. features due to the inclusion of new variables. 

In this contribution we propose use of Granger causality to analyze systems of coupled oscillators. To this aim, we 
adapt the approach of fllj ] to handle circular variables; we show, by means of several examples, that the interpretation 
of Granger causality, in terms of flow of information, is sound also in this case. 
, Let us consider a system of n interacting circular variables #i(t), 6*2 (i) , ■ ■ •, 9 n (t) discretely sampled in time. First 
^ ' of all we specify the input and output variables of the regression model. Inspired by previous works concentrating on 
t-H \ detecting the direction of coupling in interacting oscillators Q, we proceed as follows. 

For t = 1, . . . , N, N being the number of samples, we call A#j(t) the phase increments 0i(t) — 8i(t — 1), and we 
consider the problem of predicting A9i on the basis of {9j(t — 5)}j=i,..., n ;<5=i,...,m) i-e. all the phase values with lags 
from one to m. A causal relationship j — > i corresponds to an improvement of the prediction, due to the knowledge of 
Oj(t— 1),. . .,9j(t — m). According to the theory of Fourier series, predictions can be made performing linear regression 
, in the feature space spanned by the variables 
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S iS {t) = sinlhBxit -Si) + --- + £ n 9 n (t - 6 n )], (2) 

for all integers £ — (£1, . . . , £ n ) S Z" and delays 5 e {1, . . . , m} n (only terms with \£i\ -\ h \£ n \ < L are actually taken 

into account, L depending on the amount of data at disposal). All the A^-dimensional vectors {C,S} and {A9} can 
be assumed to have zero mean and unit norm, after suitable linear transformations. Now we consider the evaluation 
of the causality C(j — > i). We denote H C $t N the linear span of all vectors {C^g, Sg g} and Hq C $l N the linear span 

of those vectors {Cpg, Sjrg} with £j = 0. Decomposing H = H$ © H , we denote {u a } the orthonormal basis of H 1 - 
fl2| and calculate 

a.' 

the sum above being over significative projections [l3j . C(j — ■> i) probes the flow of information from 9j to 9i. 



C a (t) = coslWt -8 1 ) + --- + £ n 9 n {t - £„)], (1) 

and 
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In the following we will show the application of the proposed analysis to data arising from the Euler discretization 
of noisy Kuramoto's equations: 



where Sji is equal to one (zero) if there is (not) a coupling from Oj to 9i. Therefore, we take L = 2 and m = 1. 

As it has been pointed out in [lij , a fundamental property of causality estimators is the ability to discern whether 
the influence between two channels is direct or mediated. Therefore, as the first example, we consider a system of 
three coupled oscillators in which the first oscillator is coupled to the second and the second to the third (no coupling 
between the first and the third oscillators). Although there is a relevant phase correlation between the first and the 
third variables, 



our approach correctly reveals that the influence between the two time series is actually mediated by the second 
oscillator (see figure 1). 

As a second example, we consider the transient phase dynamics of a fully connected system of three oscillators. 
We fix the values of couplings so that, in the absence of noise, full synchronization of the three phases would arise 
in the large time limit. At t = 1 the phases are randomly assigned by uniform sampling in [0, 2ir). Subsequently the 
system undergoes a transient process during which the phases organize and the cross-trial phase correlation between 
oscillators, R lj (t) = | ( e ^(f ^ C*)-^ (*)) } | where the average is over initial phases, increases (figure 2-top); the asymptotic 
value of R is not exactly one due to the presence of noise. In figure 2-bottom the same process is described in terms 
of the causality between oscillators (note that in this case the different samples, needed to evaluate causality at fixed 
time, arise from many realizations of the process with varying initial phases). We observe that the causality is larger 
at beginning, when oscillators are organizing. The value of the causality at large times is much smaller. In the 
absence of noise causalities would vanish, at large t: here noise frequently perturbs the system and drives it out of 
the synchronized state. A similar cross-trial analysis, in terms of causalities, may also be performed to study the 
transient after a stimulus [l5j ]. 

We consider now the case of a binary matrix sji corresponding to an undirected graph, made of 50 nodes and 50 
links, built by means of the preferential attachment procedure [161 ]: typically, in these networks there are few nodes 
with a large number of connections (hubs), whilst most of the nodes have small connectivity. We assume that all 
oscillators in the network have the same natural frequency u>i = to: we find that the average phase synchronization 
between one node and one of its neighbors is almost independent of the connectivity k (figure 3-top). In figure 
3-middle we depict the sum of the outgoing causalities from a node, as a function of the number of its links k: it 
tends to saturate. In figure 3-bottom we depict the sum of the incoming causalities of a node as a function of the 
connectivity k: using the interpretation of causality in terms of information, these plots show that nodes with k > 4 
receive more information than they transmit and suggest the presence of a maximum amount of information that a 
node can transmit in this system. This result can be seen from the point of view of the law of diminishing marginal 
returns [13], which states that when the amount of a variable resource is increased, while other resources are kept 
fixed, the resulting change in the output will eventually diminish. In the language of social networks, this is due to 
the fact that there is a limit to the information a person may handle Turning to consider an undirected graph 
with 50 nodes, each connected to 4 randomly chosen nodes, we assign different natural frequency at each node: as 
displayed in figure 4, we find that both the average phase synchronization with neighbors and the average causality 
with neighbors are nearly independent of cjj. 

As a further application, we consider now a model of interacting thalamocortical neuronal ensembles proposed in 
[T^ | to account for the behaviour of 5 and 9 waves during anaesthesia. The model consists of three ensembles of 
Kuramoto's oscillators, namely the cortical (CO), the thalamocortical relay neurons (TC) and the thalamic reticular 
neurons (RE), having different mean natural frequencies and characterized by intra-ensemble and inter-ensemble 
coupling parameters; CO neurons receive sensory inputs from TC neurons, TC neurons receive sensory inputs from 
both CO and RE neurons, whilst RE neurons receive sensory inputs only from TC neurons. In order to simulate the 
transition from the deeply to the lightly anaesthetized state, the intra and inter couplings are swept linearly to mimic 
the effect of decreasing concentration of anaesthetic agent. In [l9[ this model has been studied by analyzing the mean 
frequencies of ensembles and the phase correlations between ensembles. In figure 5 we depict the phase correlations 
and the causalities among the three groups of neurons, as a function of /3, the overall factor of couplings. Figure 
5-top is in agreement with the results reported in [l9j , whilst figures 5-middle and 5-bottom show that the transition 
from deep to light anaesthesia is characterized by asymmetries in the causality relationships: CO neurons drive TC 
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neurons, and TC neurons drive RE neurons. These findings complete the analysis reported in [l9j . whose aim was to 
tackle the problem of anaesthetic awareness. 

Summarizing, we have proposed use of Granger causality for the analysis of systems of circular variables. Processing 
data from simulated Kuramoto systems, we have shown that causality (i) discerns direct and mediated interactions 
and (ii) is suitable to study transient phenomena. Our results, on systems living on networks built by preferential 
attachment, on one side support the interpretation of Granger causality in terms of flow of information, on the other 
hand they show that in these systems there is a maximum amount of information that an oscillator can handle, in 
accordance with the law of diminishing marginal returns. Finally, we have analyzed a recently proposed model of 
the transition from deep to light anaesthesia, and used causality to put in evidence the drive-response relationships 
between ensembles. We believe that the proposed approach will be useful to analyze real data in the form of phases. 
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FIG. 1: Causalities in a system of three coupled oscillators, with natural frequencies o»i = 0.3, ll>2 = 0.4, u>z = 0.5. The noise 
term has autocorrelation {£i(t)£j(t')} = 2DSijS(t — t'), with D = 0.01. Couplings with strength f3 are introduced from oscillator 
1 to oscillator 2 and from 2 to 3. Our estimate of the causalities are performed on time series of length N = 200, recorded in 
the steady state after the transient. The interaction 1 — » 3 is correctly recognized as mediated by oscillator 2. These results 
are robust to changes in parameters of the model. 
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FIG. 2: Phase correlation and causality, as a function of time, for a system of three oscillators, with natural frequencies 
loi = 0.25, u>2 = 0.2, UI3 = 0.15, noise strength D — 0.01 and all-to-all couplings of strength (3 = 0.3. Different samples, at fixed 
t, arise from 500 different runs of the model with random initial phases. The values of R and C, here plotted, are averaged 
over all pairs of oscillators. 
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FIG. 3: Kuramoto systems living on undirected graphs built by preferential attachment, are analyzed in terms of causality. 
Graphs have 50 nodes and 50 links, all oscillators have the same natural frequency u) = 0.4, couplings (3 = 0.3 are set between 
oscillators connected by a link. The noise term is D = 0.1. Quantities refer to the steady state, after the transient. Results are 
averaged over 20000 different networks. (Top) The mean phase correlation between a node and its neighbors is displayed as a 
function of k, the number of links of that node. It appears to be independent of k (Middle) For a node of connectivity k, the 
sum of outgoing causalities from that node to its neighbors, C ou t, is depicted. It saturates at large k. (Bottom) For a node of 
connectivity k, the sum of incoming causalities from its neighbors to that node, C; n , is depicted. For k < 4 (k > 4),d n < Cout 
(Cin > Cout)- These results are robust to changes in parameters of the model. 
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FIG. 4: Phase correlation (top) and causality (bottom) are plotted as a function of the natural frequency u) for a Kuramoto 
system living a random graph of 50 nodes, each node being connected to 4 other random nodes. Coupling is f3 — 0.05 and 
D = 0.1. Quantities refer to the steady state, after the transient. Results are averaged over 10000 different networks. Both R 
and C appear to be almost independent of to. 



6 




FIG. 5: A model to describe the transition from deep to light anaesthesia is simulated. It consists of three groups of neurons, 
CO-TC-RE, each made of 6 oscillators. The parameters (natural frequencies, inter and intra couplings) are fixed as in [l9l ]. 
As a function of an overall factor f3 (multiplying all couplings and thus simulating the effect of decreasing concentration of 
anaesthetic agent), the phase correlation between groups is depicted (top), the causalities between CO and TC groups (middle) 
and causalities between TC and RE groups (bottom). The causality from group A to group B is estimated as follows. For each 
oscillator in B, we estimate the improvement in prediction due to the inclusion, in the regression model, of all the oscillators in 
A; then we average the outcome over all the oscillators in B. 



